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Abstract. Double-mode Cepheid behaviour may arise either from non-resonant or resonant mode interaction. Inclusion of 
turbulent convection into pulsation codes by Kollath et al. led to stable double-mode F+10 Cepheid pulsation [1]. However, 
our recent computations indicate, that these models resulted from incorrect neglect of negative buoyancy effects [2]. Once 
these effects are taken into account, only some resonant double-mode models could be found. These models do not solve 
the puzzle of double-mode phenomenon, as they are restricted to narrow parameter ranges. For majority of the observed 
double-mode Cepheids non-resonant mechanism has to be operational. 
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INTRODUCTION 

Most of classical Cepheids are single-mode pulsators, 
pulsating predominantly either in the fundamental mode 
(F) or in the first overtone (lO). Few hundreds of double- 
mode variables are also known. Most of these pulsate 
simultaneously in two consecutive radial modes, either in 
fundamental and first overtone (F+IO) or in two lowest 
order overtones (10+20). 

The mechanism responsible for the double-mode be- 
haviour is far from being understood. Necessary condi- 
tion for double-mode pulsation to occur is simultaneous 
linear instability of two pulsation modes of interest. We 
consider fundamental and first overtone modes assum- 
ing for simplicity that all other modes are linearly stable. 
Fundamental mode is linearly unstable (its linear growth 
rate, 70, is positive) at the red side of the instability strip 
(IS), while first overtone is linearly unstable (71 > 0) at 
the blue side of the IS. In the middle of the IS lays the 
potential domain of the double-mode pulsators. As the 
single-mode Cepheid enters this domain during its evolu- 
tion, its pulsation state (mode selection) is determined by 
nonlinear effects. The important quantity to follow is the 
stability of the single-mode full amplitude (limit cycle) 
pulsation, with respect to perturbation in the other lin- 
early excited mode. This is expressed through the single- 
mode stability coefficients, 7 ^, which describe the sta- 
bility of the limit cycle pulsation in mode j with respect 
to perturbation in mode i. Positive value of the stability 
coefficient means that respective single-mode pulsation 
is unstable and will tend to switch into the other mode. 
Simultaneous instability of two limit cycles, say, 71,0 > 
and 70,1 > 0, leads unavoidably to F+IO double-mode 
pulsation. 



There are two mechanism that can lead to simulta- 
neous instability of both limit cycles, non-resonant and 
resonant [3]. In non-resonant mechanism two pulsation 
modes compete with each other to saturate the pulsation 
instability. If none of the modes is able to saturate the in- 
stability alone, double-mode pulsation arises. In resonant 
mechanism one of the modes (or both) is resonantly cou- 
pled with linearly damped mode. The damped, parasite 
mode, can limit the amplitude of the otherwise dominant 
mode and thus, decrease the stability of its limit cycle. If 
the mode is destabilized in a parameter range in which 
the limit cycle of the other linearly excited mode is also 
unstable, double-mode pulsation arises. 

Low order resonances among pulsation modes are rare 
in a period range occupied by double-mode variables 
and can affect the mode selection only in a rather nar- 
row parameter range. Therefore, in most of the observed 
double-mode Cepheids (both F+IO and 10+20) non- 
resonant mechanism should be operational. For many 
years, nonlinear radiative computations failed to repro- 
duce the double-mode Cepheid behaviour. It was the 
inclusion of turbulent convection into the model equa- 
tions that led to success [1]. However, our recent com- 
putations [2] indicate, that the computed non-resonant 
double-mode models resulted from incorrect physical as- 
sumption adopted by Kollath et al. [1]. We discuss these 
results below, first, providing a brief overview of turbu- 
lent convection recipes used in pulsation hydrocodes. 



TURBULENT CONVECTION RECIPES 

There are several models of turbulent convection and its 
interaction with pulsation. However, only the simplest 
models are suitable for nonlinear computations. In these 



models, generation of turbulent energy, e t , is described 
by one additional equation, which has to be solved to- 
gether with momentum and energy equations. In the 
KuhfuB model [4], which we will discuss here, these 
equations are following (see [5]): 
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Velocity u is time derivative of radius, u = dr/dt, M r is 
mass enclosed in radius r, V is specific volume, p is gas 
pressure, p t is turbulent pressure, Lj, L c and L t are radia- 
tive, convective and turbulent luminosities respectively. 
U q and E q are viscous momentum and energy transfer 
rates, respectively, and C is the term coupling the turbu- 
lent and internal energy equations. Its form is following: 



C = S D D r . 



(4) 



S is the turbulent source function, D is turbulent dis- 
sipation function and D r is radiative cooling function. 
The model is a simple phenomenological recipe, con- 
taining 8 order of unity free parameters. Nevertheless, 
all the terms have clear physical interpretation. Turbu- 
lent source function represents the buoyant forces act- 
ing on the convective eddy, responsible for its accelera- 
tion. It is proportional to superadiabatic gradient, S ~ Y, 
Y = V — V a . According to Schwarzschild criterion for 
convective instability, in stellar layers in which Y > 0, 
buoyancy represented in the source term drives the con- 
vective motions. In case of convective stability, Y < 0, 
turbulent source function is negative and contributes to 
the damping of turbulent energies. Turbulent dissipation 
term models the turbulent cascade - the decay of tur- 
bulent eddies into smaller and smaller scales down to 
molecular scale where energy is dissipated. Radiative 
cooling terms damps the turbulent energies also, as it rep- 
resents the radiative losses during the motion of the eddy. 
Viscous momentum and energy transfer rates describe 
the coupling between turbulent motions and mean, pulsa- 
tion induced gas motion. Turbulent energy is always gen- 
erated at the cost of pulsation energy, thus, eddy-viscous 
damping is an important factor limiting the pulsation am- 
plitude. For more detailed description of the model, we 
refer the reader to [5]. 

Described model was adopted in Warsaw pulsation 
codes [5]. Following notation introduced by [5] we de- 
note this model, default in all our computations, by NN. 
Very similar model was adopted in Florida-Budapest hy- 
drocode (eg. [1], [6]), which was used to compute the 
first robust double-mode Cepheid models. In the model 




FIGURE 1. Hydrodynamic integrations for particular 
Cepheid model with different initial conditions. 



adopted in that code, superadiabatic gradient, Y, was re- 
stricted to positive values. Consequently turbulent source 
function, S ~ Y+, was set equal to zero in convectively 
stable regions. Such assumption accounts to the neglect 
of a very important physical effect acting in convectively 
stable regions, namely to the neglect of negative buoy- 
ancy. In our opinion such neglect is not justified. Buoy- 
ant forces do not vanish in convectively stable regions. 
In fact, it is the negative buoyancy that restores the con- 
vective stability. Discussed assumption leads to different 
convection model (PP in the following), which we im- 
plemented in our codes as an option. Detailed compari- 
son of models computed with NN and PP convection re- 
vealed, that the double-mode behaviour observed in PP 
models resulted entirely from unphysical neglect of neg- 
ative buoyancy effects. 



NON-RESONANT DOUBLE-MODE 
CEPHEID MODELS - REVISITED 

To search and analyze the mode selection we adopted 
the methods described by [6] (see Fig. 1). For several 
models along a sequence of eg. constant luminosity and 
varying temperature, direct nonlinear integrations were 
performed. Each model was kicked with several initial 
conditions and time evolution of mode amplitudes was 
followed through the analytical signal method. Resulting 
trajectories for particular model adopting NN convection 
are shown with dashed lines in Fig. 1. Model trajectories 
were fitted with fifth order amplitude equations, which 
allowed to find all the possible pulsation states and their 
stability. Stable pulsation states (attractors) are marked 
with solid squares in Fig. 1, while unstable with open 
squares. Using interpolation, modal selection may be 
found at each temperature along a sequence. For detailed 
description of the method, we refer the reader to [2]. 



In Fig. 1 two sets of trajectories are plotted. Dashed 
lines correspond to model computed with NN convec- 
tion, while solid lines, to model of exactly the same 
parameters, but computed with PP convection. Most 
striking difference concerns the mode selection. For PP 
model, stable double-mode solution exists, while for NN 
model of the same parameters only fundamental mode 
limit cycle is stable. The other difference is a significant 
reduction of mode amplitudes of PP model, as compared 
to the NN model. Particularly, amplitude of the funda- 
mental mode is reduced very strong. Described differ- 
ences are typical. They do not depend on exact param- 
eters of the model, either convective or physical. Partic- 
ularly, along model sequences computed using PP con- 
vection, stable double-mode solutions are easily found, 
while no such solution was found for models computed 
using NN convection, despite extensive search [2]. As 
analyzed by [2] and briefly summarized below, the dif- 
ferences between NN and PP models result from the ne- 
glect of negative buoyancy effects in PP model. As such 
neglect is not justified, stable double-mode behaviour ob- 
served in PP models is not correct. 

We focus our attention on the simplest convective 
model, that is we neglect the turbulent flux, turbulent 
pressure and radiative losses. We stress however, that 
inclusion of these effects does not change the qualitative 
conclusions presented below. With these simplifications, 
turbulent energy equation is following: 

de, 

-^=S-D + E q (5) 

The crucial difference between NN and PP models ap- 
pears in convectively stable regions, where Y < 0. In 
these regions of NN model, turbulent source function is 
negative and damps the turbulent energies. In PP model 
it is set equal to zero. In both models, turbulent dissipa- 
tion term, D, always damps the turbulent energies, while 
eddy-viscous term, E q , always drives the turbulent en- 
ergies. Consequences of different treatment of turbulent 
source function are presented in Figs. 2 and 3. 

In Fig. 2 we plot the profiles of turbulent energy dur- 
ing one cycle of fundamental mode pulsation. Turbulent 
energies are plotted in logarithmic scale to highlight the 
differences in convectively stable regions, barely visible 
in a linear scale. Convectively stable, internal layers ex- 
tend below zone 70 in both NN and PP models. Here 
the differences in turbulent energy profiles are visible. 
In NN model, turbulent energies fall rapidly as supera- 
diabatic gradient changes its sign. Turbulent motions are 
braked effectively by negative buoyancy represented in 
the source term. In internal layers of NN model, tur- 
bulent energies are negligible (but nonzero, see [2]). In 
case of PP model situation is different. As superadia- 
batic gradient changes its sign, turbulent energies fall, 
but only by three orders of magnitude, as compared to 



the values in the center of the convection zone. Then, at 
the absence of negative source function, the balance be- 
tween the turbulent dissipation and eddy-viscous forces 
sets the turbulent energies at a level of 10 8 — 10 10 erg/g, 
in a large region, extending to more than 6 local pres- 
sure scale heights below the envelope convective zone 
(zones 40-70). These turbulent energies are generated 
at the cost of pulsation through the eddy-viscous forces. 
Although these energies seem small, their effect on pul- 
sation dynamics is significant, which is best visible in 
work integrals displayed in Fig. 3. In the internal layers 
of the PP model, strong eddy-viscous damping is present 
(arrows in Fig. 3). No such damping is visible in con- 
vectively stable layers of NN models, as turbulent ener- 
gies are not significant there. As a result amplitudes of 
pulsation modes computed with PP convection are much 
smaller, which is clearly visible in Fig. 1. Amplitude of 
the fundamental mode is always reduced stronger, than 
the amplitude of the first overtone due to the different 
properties of these modes in the interior of the model 1 . 

Strong reduction of fundamental mode amplitude has 
significant effect on its limit cycle stability (described by 
7i o). Lower the amplitude of the mode, more unstable 
it is against perturbation in the other mode. The run of 
stability coefficients of the fundamental and first over- 
tone limit cycles, along a sequence of models of con- 
stant luminosity is shown if Fig. 4. The qualitative sce- 
nario for both NN and PP models is similar and con- 
sistent with what we expect from non-resonant mode 
interaction. First overtone becomes more unstable (yo.i 
increases) toward the red, while fundamental mode be- 
comes more unstable (y^o increases) toward the blue. 
For NN model sequence fundamental mode limit cycle 
becomes stable at high temperatures, at which first over- 
tone limit cycle is still firmly stable. No stable-double 
mode solution is possible. Instead, we observe a mode 
selection typical for radiative computations: domains of 
single-mode first overtone pulsation (blue side of the IS) 
and fundamental mode pulsation (red side) are separated 
by either-or domain in which pulsation in either mode 
is possible, depending on initial conditions (direction of 
evolution). Situation is very different for PP model se- 
quence. Due to its strongly reduced amplitude, funda- 
mental mode is more unstable in PP sequence. It remains 
unstable at relatively low temperatures, at which first 
overtone is already unstable. Consequently in a particular 
temperature range both limit cycles are unstable, giving 
rise to stable double-mode pulsation. However, the cause 



1 Below the pulsation node in the velocity profile of the first overtone 
mode, eddy-viscous driving is very small, as E q term is proportional 
to the square of spatial derivative of velocity. Consequently, internal 
eddy-viscous damping in convectively stable regions of the PP model 
is weaker for the first overtone mode (and overtones in general). 




FIGURE 2. Profiles of turbulent energy during one cycle of 
full amplitude fundamental mode pulsation. 

of the observed double-mode behaviour is unphysical. It 
arose due to unphysical neglect of negative buoyancy ef- 
fects. 



RESONANT DOUBLE-MODE CEPHEID 
MODELS 

Despite rather disappointing results presented in the pre- 
ceding Section, we have found several stable double- 
mode models with NN convection, including negative 
buoyancy effects. In these models resonant mechanism 
is operational. Computed models include both F+IO 
Cepheids, with metallicities corresponding to Galaxy, 
and 10+20 Cepheids found during the computations in- 
tended to model the double-overtone Cepheids in the 
LMC. In the former case (F+IO), the three-mode res- 
onance of the type 2(0i = coo + 02 is the cause of the 
double-mode pulsation, while in the latter case (10+20), 
the 2:1 resonance, 2o)\ = 0)5 is operational. Details of the 
models will be published in the nearest future (Smolec & 
Moskalik, in prep.), here we present the preliminary re- 
sults. 

In Fig. 5 we display the results of numerical inte- 
gration of one Cepheid model lying close to the three- 
mode resonance center. Stable double-mode attractor is 
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FIGURE 3. Nonlinear work integrals for full amplitude fun- 
damental mode pulsation. 
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FIGURE 4. The run of stability coefficients along a sequence 
of Cepheid models of constant luminosity. 

clearly visible. It coexists with stable fundamental mode 
attractor and this mode selection is a rule in all double- 
mode models of this type that we found. The model is 
very close to the resonance center, the proximity param- 
eter, A = 2c0i/(c0o + C02) = 0.9995. Period and period 
ratios for this model are given in Fig. 5 and they are 
also typical. Somewhat too low period ratio, as com- 
pared with observation, can be compensated by slight de- 
crease of model metallicity. In all our model sequences 
double-mode behaviour was restricted to a narrow tem- 
perature range (at most 50 K, usually below 20 K). Com- 
puted double-mode models always fulfilled two condi- 
tions: (i) they were always located very close to the res- 
onance center and (z'z) in between the single-mode fun- 
damental mode pulsation domain (to the red) and F/IO 
either-or domain (to the blue), which is rather typical for 
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FIGURE 5. Hydrodynamic integrations for particular 
Cepheid model crossing the 2a>i = a>o + (Oi resonance center. 



FIGURE 6. Petersen diagram for the double-overtone LMC 
Cepheids. Different symbols indicate the mode selection in- 
ferred from nonlinear computations. 



non-resonant models. At the moment, we do not have ex- 
planation why (ii) is necessary for the double-mode be- 
haviour to occur. 

Other example of the double-mode behaviour was 
found during the model survey intended to reproduce 
the observed double-overtone Cepheids in the LMC. Our 
results are summarized in Fig. 6, showing the Petersen 
diagram for the double-overtone LMC Cepheids [7]. 
Mode selection was computed along evolutionary tracks 
of 2.5M S , 3.0M S and 3.5M S , corresponding to the first 
crossing of the IS (three leftmost tracks in Fig. 6) and 
along two horizontal paths for models with 3.0Ms and 
3.5Ms and luminosity increased by AlogL = 0.4 relative 
to evolutionary track (remaining two tracks in Fig. 6). 
Such luminosity increase was intended to model the core 
helium burning phase (blue loop). For detailed results 
of linear modeling and challenges it represent for stellar 
evolution theory, we refer the reader to [8], here focusing 
on mode selection only. 

It is clearly visible that inferred modal selection 
disagree with observation. Double-mode models were 
found only at the shortest periods. Linear analysis indi- 
cate, that the origin of the double-mode pulsation is most 
likely connected with the 2:1 resonance between the ex- 
cited first overtone and linearly damped fifth overtone. 
Indeed, resonance center, indicated with arrow in Fig. 6 
lays in the middle of the double-mode domain. Resonant 
mechanism can be operational only at shorter periods. At 
longer periods (logPi > —0.1) no low order resonances 
between pulsation modes are found and a non-resonant 
mechanism has to be operational. As in case of F+IO 
variables, no non-resonant 10+20 models were found. 

CONCLUSIONS 



works of the Florida-Budapest group (eg. [1], [6], Buch- 
ler, these proceedings). We have shown that their non- 
resonant F+IO double-mode Cepheid models resulted 
from the incorrect neglect of negative buoyancy effects. 
When these effects were included in model equations, 
no non-resonant double-mode behaviour was found. Al- 
though we have found some resonant double-mode mod- 
els including negative buoyancy in model computations, 
these models do not offer the solution to the problem 
of modelling the double-mode Cepheids. For majority 
of the stars a non-resonant mechanism has to be op- 
erational. It seems that the better treatment of convec- 
tion/pulsation coupling is necessary to solve the puzzle 
of double-mode Cepheids. 
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Presented results are rather disappointing and signif- 
icantly change the picture emerging from the recent 



